library(tidyverse)
library(plotly)
library(sf)
library(mapview)
library(tigris)
library(censusapi)
library(leaflet)
library(lehdr)
library(usmap)
library(lmtest)
library(pracma)
library(lmtest)
library(forecast)
library(vars)
library(rvest)
library(RSelenium)
library(seleniumPipes)
library(dLagM)
options(
tigris_class = "sf",
tigris_use_cache = TRUE
)
Get the case data, first for SCC.
# remDr <- remoteDriver(
# remoteServerAddr = "192.168.86.25",
# port = 4445L
# )
# remDr$open()
#
# remDr$navigate("https://app.powerbigov.us/view?r=eyJrIjoiZTg2MTlhMWQtZWE5OC00ZDI3LWE4NjAtMTU3YWYwZDRlOTNmIiwidCI6IjBhYzMyMDJmLWMzZTktNGY1Ni04MzBkLTAxN2QwOWQxNmIzZiJ9")
#
# webElem <- remDr$findElements(using = "class", value = "column")
#
# cases <-
# 1:length(webElem) %>%
# map(function(x){
# webElem[[x]]$getElementAttribute("aria-label") %>% as.character()
# }) %>%
# unlist() %>%
# as.data.frame()
#
# scc_cumulative_cases <-
# cases %>%
# rename(text = ".") %>%
# filter(grepl("Total_cases",text)) %>%
# separate(text, c("date","cases"), sep = "\\.") %>%
# mutate(
# date =
# substr(date,6,nchar(.)) %>%
# as.Date("%A, %B %d, %Y"),
# cases =
# substr(cases,13,nchar(.)) %>%
# as.numeric()
# )
#
# saveRDS(scc_cumulative_cases, file = "/Users/simonespeizer/Documents/2020 Spring Quarter/CEE 218Z/covid19/Simone_Speizer/scc_cumulative_cases.rds")
scc_cumulative_cases <- readRDS("/Users/simonespeizer/Documents/2020 Spring Quarter/CEE 218Z/covid19/Simone_Speizer/scc_cumulative_cases.rds")
Also for SMC.
# remDr$navigate("https://app.powerbigov.us/view?r=eyJrIjoiMWI5NmE5M2ItOTUwMC00NGNmLWEzY2UtOTQyODA1YjQ1NWNlIiwidCI6IjBkZmFmNjM1LWEwNGQtNDhjYy1hN2UzLTZkYTFhZjA4ODNmOSJ9")
#
# webElem <- remDr$findElements(using = "class", value = "column")
#
# tests <-
# 1:length(webElem) %>%
# map(function(x){
# webElem[[x]]$getElementAttribute("aria-label") %>% as.character()
# }) %>%
# unlist() %>%
# as.data.frame()
#
# tests_clean <-
# tests %>%
# rename(text = ".") %>%
# separate(text, c("date","test_text"), sep = "\\.") %>%
# separate(test_text, c(NA,"type",NA,"tests")) %>%
# mutate(
# date =
# substr(date,23,nchar(.)) %>%
# as.Date("%A, %B %d, %Y"),
# tests =
# tests %>%
# as.numeric()
# ) %>%
# spread(
# key = type,
# value = tests
# ) %>%
# mutate(
# total = Negative + Pending + Positive,
# perc_positive = Positive/total,
# perc_positive_movavg = movavg(perc_positive, 7, type = "s")
# )
#
# smc_cumulative_cases <- tests_clean %>%
# mutate(cumulative_cases = cumsum(Positive), cumulative_negative = cumsum(Negative), cumulative_total = cumulative_cases+cumulative_negative, perc_positive_cumulative = cumulative_cases*100 / cumulative_total, perc_positive_cumulative_mov = movavg(perc_positive_cumulative, 7, type = "s"))
#
# saveRDS(smc_cumulative_cases, file = "/Users/simonespeizer/Documents/2020 Spring Quarter/CEE 218Z/covid19/Simone_Speizer/smc_cumulative_cases.rds")
smc_cumulative_cases <- readRDS("/Users/simonespeizer/Documents/2020 Spring Quarter/CEE 218Z/covid19/Simone_Speizer/smc_cumulative_cases.rds")
Get social distancing data.
scc_blockgroups <-
block_groups("CA","Santa Clara", cb=F, progress_bar=F) %>%
st_transform('+proj=longlat +datum=WGS84')
smc_blockgroups <-
block_groups("CA","San Mateo", cb=F, progress_bar=F) %>%
st_transform('+proj=longlat +datum=WGS84')
bay_sd <- readRDS("/Users/simonespeizer/pCloud Drive/Shared/SFBI/Restricted Data Library/Safegraph/covid19analysis/bay_socialdistancing_v2.rds") %>%
mutate(date = date_range_start %>% substr(1,10) %>% as.Date())
# obtaining weekends
weekends <- bay_sd %>%
filter(!duplicated(date)) %>%
arrange(date) %>%
mutate(weekend = ifelse((date %>% as.numeric()) %% 7 %in% c(2,3), T, F)) %>%
dplyr::select(date,weekend)
bay_sd <- bay_sd %>% left_join(weekends)
scc_cases_sd_daily <- bay_sd %>%
filter(origin_census_block_group %in% scc_blockgroups$GEOID) %>%
group_by(date) %>%
summarize(total_at_home = sum(completely_home_device_count), total_devices = sum(device_count)) %>%
mutate(
percent_at_home = total_at_home*100/total_devices,
percent_leaving_home = (100 - percent_at_home),
) %>%
left_join(
scc_cumulative_cases
) %>%
filter(date >= min(scc_cumulative_cases$date))
# get the baseline percent of people at home
pre_case_growth_at_home_scc <- bay_sd %>%
filter(date < min(scc_cumulative_cases$date)) %>%
filter(origin_census_block_group %in% scc_blockgroups$GEOID) %>%
summarize(total_at_home = sum(completely_home_device_count), total_devices = sum(device_count)) %>%
mutate(percent_at_home = total_at_home*100/total_devices, percent_leaving_home = (100 - percent_at_home))
# include change in percent leaving home
scc_cases_sd_daily <- scc_cases_sd_daily %>%
mutate(leaving_home_dif = percent_leaving_home - pre_case_growth_at_home_scc$percent_leaving_home[1],
log_cases = log(cases))
# compute number of differences for stationarity
ndiffs(scc_cases_sd_daily$cases)
## [1] 2
ndiffs(scc_cases_sd_daily$log_cases[-1])
## [1] 2
ndiffs(scc_cases_sd_daily$leaving_home_dif)
## [1] 1
scc_test_big <- scc_cases_sd_daily %>%
mutate(
dlog_cases = c(NA,diff(log_cases)),
d2log_cases = c(NA,NA,diff(log_cases, differences = 2)),
dcases = c(NA,diff(cases)),
d2cases = c(NA,NA,diff(dcases, differences = 2)),
dleaving = c(NA,diff(leaving_home_dif)),
d2leaving = c(NA,NA,diff(leaving_home_dif, differences = 2)),
cases_mov = movavg(cases, 7, type = "s"),
log_cases_mov = movavg(log_cases, 7, type = "s"),
dlog_cases_mov = c(NA,diff(log_cases_mov)),
d2log_cases_mov = c(NA,NA,diff(log_cases_mov, differences = 2)),
dcases_mov = c(NA,diff(cases_mov)),
d2cases_mov = c(NA,diff(dcases_mov)),
leaving_mov = movavg(leaving_home_dif, 7, type = "s"),
dleaving_mov = c(NA,diff(leaving_mov)),
d2leaving_mov = c(NA,diff(dleaving_mov)),
leaving4 = c(rep(NA,4), scc_cases_sd_daily$leaving_home_dif[1:(nrow(scc_cases_sd_daily)-4)]),
dleaving4 = c(NA,diff(leaving4)),
d2leaving4 = c(NA,NA,diff(leaving4, differences = 2)),
leaving3 = c(rep(NA,3), scc_cases_sd_daily$leaving_home_dif[1:(nrow(scc_cases_sd_daily)-3)]),
leaving3_mov = movavg(leaving3, 7, type = "s"),
dleaving3_mov = c(NA,diff(leaving3_mov)),
d2leaving3_mov = c(NA,NA,diff(leaving3_mov, differences = 2)),
leaving4_mov = movavg(leaving4, 7, type = "s"),
dleaving4_mov = c(NA,diff(leaving4_mov)),
d2leaving4_mov = c(NA,NA,diff(leaving4_mov, differences = 2)),
leaving5 = c(rep(NA,5), scc_cases_sd_daily$leaving_home_dif[1:(nrow(scc_cases_sd_daily)-5)]),
leaving5_mov = movavg(leaving5, 7, type = "s"),
dleaving5_mov = c(NA,diff(leaving5_mov)),
d2leaving5_mov = c(NA,NA,diff(leaving5_mov, differences = 2)),
leaving6 = c(rep(NA,6), scc_cases_sd_daily$leaving_home_dif[1:(nrow(scc_cases_sd_daily)-6)]),
leaving6_mov = movavg(leaving6, 7, type = "s"),
dleaving6_mov = c(NA,diff(leaving6_mov)),
d2leaving6_mov = c(NA,NA,diff(leaving6_mov, differences = 2)),
leaving7 = c(rep(NA,7), scc_cases_sd_daily$leaving_home_dif[1:(nrow(scc_cases_sd_daily)-7)]),
leaving7_mov = movavg(leaving7, 7, type = "s"),
dleaving7_mov = c(NA,diff(leaving7_mov)),
d2leaving7_mov = c(NA,NA,diff(leaving7_mov, differences = 2)),
leaving8 = c(rep(NA,8), scc_cases_sd_daily$leaving_home_dif[1:(nrow(scc_cases_sd_daily)-8)]),
leaving8_mov = movavg(leaving8, 7, type = "s"),
dleaving8_mov = c(NA,diff(leaving8_mov)),
d2leaving8_mov = c(NA,NA,diff(leaving8_mov, differences = 2)),
leaving9 = c(rep(NA,9), scc_cases_sd_daily$leaving_home_dif[1:(nrow(scc_cases_sd_daily)-9)]),
leaving9_mov = movavg(leaving9, 7, type = "s"),
dleaving9_mov = c(NA,diff(leaving9_mov)),
d2leaving9_mov = c(NA,NA,diff(leaving9_mov, differences = 2)),
leaving10 = c(rep(NA,10), scc_cases_sd_daily$leaving_home_dif[1:(nrow(scc_cases_sd_daily)-10)]),
leaving10_mov = movavg(leaving10, 7, type = "s"),
dleaving10_mov = c(NA,diff(leaving10_mov)),
d2leaving10_mov = c(NA,NA,diff(leaving10_mov, differences = 2)),
past_cases = c(NA, scc_cases_sd_daily$cases[1:(nrow(scc_cases_sd_daily)-1)]),
cases_growth_daily = (dcases / past_cases)*100,
cases_growth_daily_mov = movavg(cases_growth_daily, 7, type = "s")
) %>%
filter(date >= "2020-03-01")
ndiffs(scc_cases_sd_daily$cases_growth_daily)
## [1] 0
Granger
grangertest(
d2log_cases_mov ~ dleaving_mov,
order = 3,
data = scc_test_big
)
## Granger causality test
##
## Model 1: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:3) + Lags(dleaving_mov, 1:3)
## Model 2: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:3)
## Res.Df Df F Pr(>F)
## 1 60
## 2 63 -3 7.7236 0.0001912 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
dleaving_mov ~ d2log_cases_mov,
order = 3,
data = scc_test_big
)
## Granger causality test
##
## Model 1: dleaving_mov ~ Lags(dleaving_mov, 1:3) + Lags(d2log_cases_mov, 1:3)
## Model 2: dleaving_mov ~ Lags(dleaving_mov, 1:3)
## Res.Df Df F Pr(>F)
## 1 60
## 2 63 -3 1.4405 0.2399
grangertest(
d2log_cases_mov ~ dleaving_mov,
order = 4,
data = scc_test_big
)
## Granger causality test
##
## Model 1: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:4) + Lags(dleaving_mov, 1:4)
## Model 2: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:4)
## Res.Df Df F Pr(>F)
## 1 57
## 2 61 -4 7.6572 5.21e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
dleaving_mov ~ d2log_cases_mov,
order = 4,
data = scc_test_big
)
## Granger causality test
##
## Model 1: dleaving_mov ~ Lags(dleaving_mov, 1:4) + Lags(d2log_cases_mov, 1:4)
## Model 2: dleaving_mov ~ Lags(dleaving_mov, 1:4)
## Res.Df Df F Pr(>F)
## 1 57
## 2 61 -4 0.8209 0.5172
grangertest(
d2log_cases_mov ~ dleaving_mov,
order = 5,
data = scc_test_big
)
## Granger causality test
##
## Model 1: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:5) + Lags(dleaving_mov, 1:5)
## Model 2: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:5)
## Res.Df Df F Pr(>F)
## 1 54
## 2 59 -5 5.7937 0.0002365 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
dleaving_mov ~ d2log_cases_mov,
order = 5,
data = scc_test_big
)
## Granger causality test
##
## Model 1: dleaving_mov ~ Lags(dleaving_mov, 1:5) + Lags(d2log_cases_mov, 1:5)
## Model 2: dleaving_mov ~ Lags(dleaving_mov, 1:5)
## Res.Df Df F Pr(>F)
## 1 54
## 2 59 -5 0.7215 0.6102
grangertest(
d2log_cases_mov ~ dleaving_mov,
order = 6,
data = scc_test_big
)
## Granger causality test
##
## Model 1: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:6) + Lags(dleaving_mov, 1:6)
## Model 2: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:6)
## Res.Df Df F Pr(>F)
## 1 51
## 2 57 -6 4.7436 0.0006565 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
dleaving_mov ~ d2log_cases_mov,
order = 6,
data = scc_test_big
)
## Granger causality test
##
## Model 1: dleaving_mov ~ Lags(dleaving_mov, 1:6) + Lags(d2log_cases_mov, 1:6)
## Model 2: dleaving_mov ~ Lags(dleaving_mov, 1:6)
## Res.Df Df F Pr(>F)
## 1 51
## 2 57 -6 0.8913 0.5083
grangertest(
d2log_cases_mov ~ dleaving_mov,
order = 7,
data = scc_test_big
)
## Granger causality test
##
## Model 1: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:7) + Lags(dleaving_mov, 1:7)
## Model 2: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:7)
## Res.Df Df F Pr(>F)
## 1 48
## 2 55 -7 2.9577 0.01166 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
dleaving_mov ~ d2log_cases_mov,
order = 7,
data = scc_test_big
)
## Granger causality test
##
## Model 1: dleaving_mov ~ Lags(dleaving_mov, 1:7) + Lags(d2log_cases_mov, 1:7)
## Model 2: dleaving_mov ~ Lags(dleaving_mov, 1:7)
## Res.Df Df F Pr(>F)
## 1 48
## 2 55 -7 1.6548 0.1431
grangertest(
d2log_cases_mov ~ dleaving_mov,
order = 8,
data = scc_test_big
)
## Granger causality test
##
## Model 1: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:8) + Lags(dleaving_mov, 1:8)
## Model 2: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:8)
## Res.Df Df F Pr(>F)
## 1 45
## 2 53 -8 5.0121 0.0001761 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
dleaving_mov ~ d2log_cases_mov,
order = 8,
data = scc_test_big
)
## Granger causality test
##
## Model 1: dleaving_mov ~ Lags(dleaving_mov, 1:8) + Lags(d2log_cases_mov, 1:8)
## Model 2: dleaving_mov ~ Lags(dleaving_mov, 1:8)
## Res.Df Df F Pr(>F)
## 1 45
## 2 53 -8 1.9809 0.07091 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
d2log_cases_mov ~ dleaving_mov,
order = 9,
data = scc_test_big
)
## Granger causality test
##
## Model 1: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:9) + Lags(dleaving_mov, 1:9)
## Model 2: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:9)
## Res.Df Df F Pr(>F)
## 1 42
## 2 51 -9 3.8325 0.00133 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
dleaving_mov ~ d2log_cases_mov,
order = 9,
data = scc_test_big
)
## Granger causality test
##
## Model 1: dleaving_mov ~ Lags(dleaving_mov, 1:9) + Lags(d2log_cases_mov, 1:9)
## Model 2: dleaving_mov ~ Lags(dleaving_mov, 1:9)
## Res.Df Df F Pr(>F)
## 1 42
## 2 51 -9 1.6247 0.1394
grangertest(
d2log_cases_mov ~ dleaving_mov,
order = 10,
data = scc_test_big
)
## Granger causality test
##
## Model 1: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:10) + Lags(dleaving_mov, 1:10)
## Model 2: d2log_cases_mov ~ Lags(d2log_cases_mov, 1:10)
## Res.Df Df F Pr(>F)
## 1 39
## 2 49 -10 3.6794 0.001578 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
dleaving_mov ~ d2log_cases_mov,
order = 10,
data = scc_test_big
)
## Granger causality test
##
## Model 1: dleaving_mov ~ Lags(dleaving_mov, 1:10) + Lags(d2log_cases_mov, 1:10)
## Model 2: dleaving_mov ~ Lags(dleaving_mov, 1:10)
## Res.Df Df F Pr(>F)
## 1 39
## 2 49 -10 1.462 0.1907
Granger for not moving average
grangertest(
d2log_cases ~ dleaving,
order = 3,
data = scc_test_big
)
## Granger causality test
##
## Model 1: d2log_cases ~ Lags(d2log_cases, 1:3) + Lags(dleaving, 1:3)
## Model 2: d2log_cases ~ Lags(d2log_cases, 1:3)
## Res.Df Df F Pr(>F)
## 1 60
## 2 63 -3 4.9045 0.004095 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
dleaving ~ d2log_cases,
order = 3,
data = scc_test_big
)
## Granger causality test
##
## Model 1: dleaving ~ Lags(dleaving, 1:3) + Lags(d2log_cases, 1:3)
## Model 2: dleaving ~ Lags(dleaving, 1:3)
## Res.Df Df F Pr(>F)
## 1 60
## 2 63 -3 1.3908 0.2543
grangertest(
d2log_cases ~ dleaving,
order = 4,
data = scc_test_big
)
## Granger causality test
##
## Model 1: d2log_cases ~ Lags(d2log_cases, 1:4) + Lags(dleaving, 1:4)
## Model 2: d2log_cases ~ Lags(d2log_cases, 1:4)
## Res.Df Df F Pr(>F)
## 1 57
## 2 61 -4 7.3475 7.628e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
dleaving ~ d2log_cases,
order = 4,
data = scc_test_big
)
## Granger causality test
##
## Model 1: dleaving ~ Lags(dleaving, 1:4) + Lags(d2log_cases, 1:4)
## Model 2: dleaving ~ Lags(dleaving, 1:4)
## Res.Df Df F Pr(>F)
## 1 57
## 2 61 -4 2.3217 0.06763 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
d2log_cases ~ dleaving,
order = 5,
data = scc_test_big
)
## Granger causality test
##
## Model 1: d2log_cases ~ Lags(d2log_cases, 1:5) + Lags(dleaving, 1:5)
## Model 2: d2log_cases ~ Lags(d2log_cases, 1:5)
## Res.Df Df F Pr(>F)
## 1 54
## 2 59 -5 5.0124 0.0007638 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
dleaving ~ d2log_cases,
order = 5,
data = scc_test_big
)
## Granger causality test
##
## Model 1: dleaving ~ Lags(dleaving, 1:5) + Lags(d2log_cases, 1:5)
## Model 2: dleaving ~ Lags(dleaving, 1:5)
## Res.Df Df F Pr(>F)
## 1 54
## 2 59 -5 1.7951 0.1294
grangertest(
d2log_cases ~ dleaving,
order = 6,
data = scc_test_big
)
## Granger causality test
##
## Model 1: d2log_cases ~ Lags(d2log_cases, 1:6) + Lags(dleaving, 1:6)
## Model 2: d2log_cases ~ Lags(d2log_cases, 1:6)
## Res.Df Df F Pr(>F)
## 1 51
## 2 57 -6 3.2324 0.00908 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
dleaving ~ d2log_cases,
order = 6,
data = scc_test_big
)
## Granger causality test
##
## Model 1: dleaving ~ Lags(dleaving, 1:6) + Lags(d2log_cases, 1:6)
## Model 2: dleaving ~ Lags(dleaving, 1:6)
## Res.Df Df F Pr(>F)
## 1 51
## 2 57 -6 1.3455 0.2546
grangertest(
d2log_cases ~ dleaving,
order = 7,
data = scc_test_big
)
## Granger causality test
##
## Model 1: d2log_cases ~ Lags(d2log_cases, 1:7) + Lags(dleaving, 1:7)
## Model 2: d2log_cases ~ Lags(d2log_cases, 1:7)
## Res.Df Df F Pr(>F)
## 1 48
## 2 55 -7 3.1376 0.008239 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
dleaving ~ d2log_cases,
order = 7,
data = scc_test_big
)
## Granger causality test
##
## Model 1: dleaving ~ Lags(dleaving, 1:7) + Lags(d2log_cases, 1:7)
## Model 2: dleaving ~ Lags(dleaving, 1:7)
## Res.Df Df F Pr(>F)
## 1 48
## 2 55 -7 1.0694 0.3975
grangertest(
d2log_cases ~ dleaving,
order = 8,
data = scc_test_big
)
## Granger causality test
##
## Model 1: d2log_cases ~ Lags(d2log_cases, 1:8) + Lags(dleaving, 1:8)
## Model 2: d2log_cases ~ Lags(d2log_cases, 1:8)
## Res.Df Df F Pr(>F)
## 1 45
## 2 53 -8 3.4128 0.0038 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
dleaving ~ d2log_cases,
order = 8,
data = scc_test_big
)
## Granger causality test
##
## Model 1: dleaving ~ Lags(dleaving, 1:8) + Lags(d2log_cases, 1:8)
## Model 2: dleaving ~ Lags(dleaving, 1:8)
## Res.Df Df F Pr(>F)
## 1 45
## 2 53 -8 1.6413 0.14
grangertest(
d2log_cases ~ dleaving,
order = 9,
data = scc_test_big
)
## Granger causality test
##
## Model 1: d2log_cases ~ Lags(d2log_cases, 1:9) + Lags(dleaving, 1:9)
## Model 2: d2log_cases ~ Lags(d2log_cases, 1:9)
## Res.Df Df F Pr(>F)
## 1 42
## 2 51 -9 2.9897 0.007637 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
dleaving ~ d2log_cases,
order = 9,
data = scc_test_big
)
## Granger causality test
##
## Model 1: dleaving ~ Lags(dleaving, 1:9) + Lags(d2log_cases, 1:9)
## Model 2: dleaving ~ Lags(dleaving, 1:9)
## Res.Df Df F Pr(>F)
## 1 42
## 2 51 -9 1.5613 0.1587
grangertest(
d2log_cases ~ dleaving,
order = 10,
data = scc_test_big
)
## Granger causality test
##
## Model 1: d2log_cases ~ Lags(d2log_cases, 1:10) + Lags(dleaving, 1:10)
## Model 2: d2log_cases ~ Lags(d2log_cases, 1:10)
## Res.Df Df F Pr(>F)
## 1 39
## 2 49 -10 2.3634 0.02696 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
dleaving ~ d2log_cases,
order = 10,
data = scc_test_big
)
## Granger causality test
##
## Model 1: dleaving ~ Lags(dleaving, 1:10) + Lags(d2log_cases, 1:10)
## Model 2: dleaving ~ Lags(dleaving, 1:10)
## Res.Df Df F Pr(>F)
## 1 39
## 2 49 -10 0.8308 0.6021
Granger with case growth rate
grangertest(
cases_growth_daily ~ dleaving,
order = 2,
data = scc_test_big
)
## Granger causality test
##
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:2) + Lags(dleaving, 1:2)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:2)
## Res.Df Df F Pr(>F)
## 1 63
## 2 65 -2 0.8535 0.4308
grangertest(
dleaving ~ cases_growth_daily,
order = 2,
data = scc_test_big
)
## Granger causality test
##
## Model 1: dleaving ~ Lags(dleaving, 1:2) + Lags(cases_growth_daily, 1:2)
## Model 2: dleaving ~ Lags(dleaving, 1:2)
## Res.Df Df F Pr(>F)
## 1 63
## 2 65 -2 3.3287 0.04224 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
cases_growth_daily ~ dleaving,
order = 3,
data = scc_test_big
)
## Granger causality test
##
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:3) + Lags(dleaving, 1:3)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:3)
## Res.Df Df F Pr(>F)
## 1 60
## 2 63 -3 2.343 0.08205 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
dleaving ~ cases_growth_daily,
order = 3,
data = scc_test_big
)
## Granger causality test
##
## Model 1: dleaving ~ Lags(dleaving, 1:3) + Lags(cases_growth_daily, 1:3)
## Model 2: dleaving ~ Lags(dleaving, 1:3)
## Res.Df Df F Pr(>F)
## 1 60
## 2 63 -3 3.5172 0.02035 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
cases_growth_daily ~ dleaving,
order = 4,
data = scc_test_big
)
## Granger causality test
##
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:4) + Lags(dleaving, 1:4)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:4)
## Res.Df Df F Pr(>F)
## 1 57
## 2 61 -4 4.3608 0.003808 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
dleaving ~ cases_growth_daily,
order = 4,
data = scc_test_big
)
## Granger causality test
##
## Model 1: dleaving ~ Lags(dleaving, 1:4) + Lags(cases_growth_daily, 1:4)
## Model 2: dleaving ~ Lags(dleaving, 1:4)
## Res.Df Df F Pr(>F)
## 1 57
## 2 61 -4 2.7221 0.03821 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
cases_growth_daily ~ dleaving,
order = 5,
data = scc_test_big
)
## Granger causality test
##
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:5) + Lags(dleaving, 1:5)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:5)
## Res.Df Df F Pr(>F)
## 1 54
## 2 59 -5 3.8974 0.004343 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
dleaving ~ cases_growth_daily,
order = 5,
data = scc_test_big
)
## Granger causality test
##
## Model 1: dleaving ~ Lags(dleaving, 1:5) + Lags(cases_growth_daily, 1:5)
## Model 2: dleaving ~ Lags(dleaving, 1:5)
## Res.Df Df F Pr(>F)
## 1 54
## 2 59 -5 3.1499 0.01444 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
cases_growth_daily ~ dleaving,
order = 6,
data = scc_test_big
)
## Granger causality test
##
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:6) + Lags(dleaving, 1:6)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:6)
## Res.Df Df F Pr(>F)
## 1 51
## 2 57 -6 4.9432 0.0004699 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
dleaving ~ cases_growth_daily,
order = 6,
data = scc_test_big
)
## Granger causality test
##
## Model 1: dleaving ~ Lags(dleaving, 1:6) + Lags(cases_growth_daily, 1:6)
## Model 2: dleaving ~ Lags(dleaving, 1:6)
## Res.Df Df F Pr(>F)
## 1 51
## 2 57 -6 2.4259 0.03866 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
cases_growth_daily ~ dleaving,
order = 7,
data = scc_test_big
)
## Granger causality test
##
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:7) + Lags(dleaving, 1:7)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:7)
## Res.Df Df F Pr(>F)
## 1 48
## 2 55 -7 4.3839 0.0007929 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
dleaving ~ cases_growth_daily,
order = 7,
data = scc_test_big
)
## Granger causality test
##
## Model 1: dleaving ~ Lags(dleaving, 1:7) + Lags(cases_growth_daily, 1:7)
## Model 2: dleaving ~ Lags(dleaving, 1:7)
## Res.Df Df F Pr(>F)
## 1 48
## 2 55 -7 1.4992 0.1903
grangertest(
cases_growth_daily ~ dleaving,
order = 8,
data = scc_test_big
)
## Granger causality test
##
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:8) + Lags(dleaving, 1:8)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:8)
## Res.Df Df F Pr(>F)
## 1 45
## 2 53 -8 3.1117 0.006981 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
dleaving ~ cases_growth_daily,
order = 8,
data = scc_test_big
)
## Granger causality test
##
## Model 1: dleaving ~ Lags(dleaving, 1:8) + Lags(cases_growth_daily, 1:8)
## Model 2: dleaving ~ Lags(dleaving, 1:8)
## Res.Df Df F Pr(>F)
## 1 45
## 2 53 -8 2.9246 0.01022 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
cases_growth_daily ~ dleaving,
order = 9,
data = scc_test_big
)
## Granger causality test
##
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:9) + Lags(dleaving, 1:9)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:9)
## Res.Df Df F Pr(>F)
## 1 42
## 2 51 -9 2.14 0.04708 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
dleaving ~ cases_growth_daily,
order = 9,
data = scc_test_big
)
## Granger causality test
##
## Model 1: dleaving ~ Lags(dleaving, 1:9) + Lags(cases_growth_daily, 1:9)
## Model 2: dleaving ~ Lags(dleaving, 1:9)
## Res.Df Df F Pr(>F)
## 1 42
## 2 51 -9 2.2421 0.03783 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
cases_growth_daily ~ dleaving,
order = 10,
data = scc_test_big
)
## Granger causality test
##
## Model 1: cases_growth_daily ~ Lags(cases_growth_daily, 1:10) + Lags(dleaving, 1:10)
## Model 2: cases_growth_daily ~ Lags(cases_growth_daily, 1:10)
## Res.Df Df F Pr(>F)
## 1 39
## 2 49 -10 2.0268 0.0567 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
grangertest(
dleaving ~ cases_growth_daily,
order = 10,
data = scc_test_big
)
## Granger causality test
##
## Model 1: dleaving ~ Lags(dleaving, 1:10) + Lags(cases_growth_daily, 1:10)
## Model 2: dleaving ~ Lags(dleaving, 1:10)
## Res.Df Df F Pr(>F)
## 1 39
## 2 49 -10 3.3384 0.003229 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Also try ardlDlm
# 4th order lag on x only
testing_ardl4 = ardlDlm(x = scc_test_big$dleaving_mov, y = scc_test_big$d2log_cases_mov, p = 4, q = 4, remove = list(p = c(0,1,2,3), q=c()))
summary(testing_ardl4)
##
## Time series regression with "ts" data:
## Start = 5, End = 70
##
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0208895 -0.0015317 0.0005453 0.0019138 0.0263403
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0012166 0.0010638 -1.144 0.257296
## X.4 -0.0007294 0.0013659 -0.534 0.595297
## Y.1 -0.0499038 0.1301945 -0.383 0.702851
## Y.2 0.2326294 0.1336737 1.740 0.086937 .
## Y.3 0.4296713 0.1104830 3.889 0.000255 ***
## Y.4 0.0512616 0.1204743 0.425 0.671996
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.007618 on 60 degrees of freedom
## Multiple R-squared: 0.2415, Adjusted R-squared: 0.1782
## F-statistic: 3.82 on 5 and 60 DF, p-value: 0.004534
testing_ardl4 = ardlDlm(x = scc_test_big$dleaving, y = scc_test_big$d2log_cases, p = 4, q = 4, remove = list(p = c(0,1,2,3), q=c()))
summary(testing_ardl4)
##
## Time series regression with "ts" data:
## Start = 5, End = 70
##
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.107955 -0.014411 0.001501 0.010923 0.095726
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.006614 0.004592 -1.441 0.154912
## X.4 0.002060 0.001416 1.455 0.150881
## Y.1 -0.681884 0.126119 -5.407 1.17e-06 ***
## Y.2 -0.571356 0.138753 -4.118 0.000119 ***
## Y.3 -0.441065 0.144300 -3.057 0.003340 **
## Y.4 -0.192237 0.105759 -1.818 0.074105 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03626 on 60 degrees of freedom
## Multiple R-squared: 0.3523, Adjusted R-squared: 0.2983
## F-statistic: 6.526 on 5 and 60 DF, p-value: 6.593e-05
# 5th order lag on x only
testing_ardl5 = ardlDlm(x = scc_test_big$dleaving_mov, y = scc_test_big$d2log_cases_mov, p = 5, q = 5, remove = list(p = c(0,1,2,3,4), q=c()))
summary(testing_ardl5)
##
## Time series regression with "ts" data:
## Start = 6, End = 70
##
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.021418 -0.001004 0.001145 0.002117 0.015301
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0018571 0.0009125 -2.035 0.0464 *
## X.5 0.0012634 0.0011535 1.095 0.2779
## Y.1 0.1035728 0.1085414 0.954 0.3439
## Y.2 -0.0214885 0.1097219 -0.196 0.8454
## Y.3 0.1739958 0.1154550 1.507 0.1372
## Y.4 -0.0700718 0.1040173 -0.674 0.5032
## Y.5 0.1172137 0.1014753 1.155 0.2528
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.006404 on 58 degrees of freedom
## Multiple R-squared: 0.2163, Adjusted R-squared: 0.1353
## F-statistic: 2.669 on 6 and 58 DF, p-value: 0.0235
testing_ardl5 = ardlDlm(x = scc_test_big$dleaving, y = scc_test_big$d2log_cases, p = 5, q = 5, remove = list(p = c(0,1,2,3,4), q=c()))
summary(testing_ardl5)
##
## Time series regression with "ts" data:
## Start = 6, End = 70
##
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.102828 -0.010813 0.003215 0.012146 0.081908
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.007695 0.004079 -1.887 0.0642 .
## X.5 -0.003154 0.001248 -2.527 0.0143 *
## Y.1 -0.599547 0.111878 -5.359 1.51e-06 ***
## Y.2 -0.303334 0.133295 -2.276 0.0266 *
## Y.3 -0.319233 0.136188 -2.344 0.0225 *
## Y.4 0.183904 0.134464 1.368 0.1767
## Y.5 0.218774 0.094149 2.324 0.0237 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03142 on 58 degrees of freedom
## Multiple R-squared: 0.5057, Adjusted R-squared: 0.4546
## F-statistic: 9.891 on 6 and 58 DF, p-value: 1.691e-07
# 6th order lag on x only
testing_ardl6 = ardlDlm(x = scc_test_big$dleaving_mov, y = scc_test_big$d2log_cases_mov, p = 6, q = 6, remove = list(p = c(0,1,2,3,4,5), q=c()))
summary(testing_ardl6)
##
## Time series regression with "ts" data:
## Start = 7, End = 70
##
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.020128 -0.002107 0.001037 0.001714 0.014682
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0019052 0.0009642 -1.976 0.0531 .
## X.6 0.0013845 0.0011791 1.174 0.2453
## Y.1 0.1002644 0.1327362 0.755 0.4532
## Y.2 0.0306193 0.1104764 0.277 0.7827
## Y.3 0.2379076 0.1110240 2.143 0.0365 *
## Y.4 -0.1005564 0.1190132 -0.845 0.4018
## Y.5 0.1042246 0.1055226 0.988 0.3275
## Y.6 -0.1164147 0.1036735 -1.123 0.2663
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.006468 on 56 degrees of freedom
## Multiple R-squared: 0.2147, Adjusted R-squared: 0.1165
## F-statistic: 2.187 on 7 and 56 DF, p-value: 0.04917
testing_ardl6 = ardlDlm(x = scc_test_big$dleaving, y = scc_test_big$d2log_cases, p = 6, q = 6, remove = list(p = c(0,1,2,3,4,5), q=c()))
summary(testing_ardl6)
##
## Time series regression with "ts" data:
## Start = 7, End = 70
##
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.08285 -0.01097 0.00638 0.01140 0.07022
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0118381 0.0036709 -3.225 0.002105 **
## X.6 0.0003069 0.0011401 0.269 0.788801
## Y.1 -0.6410868 0.1137748 -5.635 5.92e-07 ***
## Y.2 -0.4546921 0.1185650 -3.835 0.000321 ***
## Y.3 -0.5892889 0.1205805 -4.887 8.96e-06 ***
## Y.4 -0.1817903 0.1235876 -1.471 0.146904
## Y.5 -0.3268229 0.1184621 -2.759 0.007821 **
## Y.6 -0.3825815 0.0853473 -4.483 3.70e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02722 on 56 degrees of freedom
## Multiple R-squared: 0.6393, Adjusted R-squared: 0.5942
## F-statistic: 14.18 on 7 and 56 DF, p-value: 1.979e-10
# 7th order lag on x only
testing_ardl7 = ardlDlm(x = scc_test_big$dleaving_mov, y = scc_test_big$d2log_cases_mov, p = 7, q = 7, remove = list(p = c(0,1,2,3,4,5,6), q=c()))
summary(testing_ardl7)
##
## Time series regression with "ts" data:
## Start = 8, End = 70
##
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.022767 -0.001574 0.000301 0.001407 0.011354
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0009224 0.0008347 -1.105 0.27400
## X.7 0.0010192 0.0009881 1.031 0.30691
## Y.1 0.0715185 0.1106610 0.646 0.52083
## Y.2 0.3304048 0.1104324 2.992 0.00417 **
## Y.3 0.1465600 0.0914905 1.602 0.11501
## Y.4 0.0871186 0.0955925 0.911 0.36616
## Y.5 0.1287338 0.0991496 1.298 0.19967
## Y.6 -0.0561249 0.0880841 -0.637 0.52670
## Y.7 -0.2286607 0.0867580 -2.636 0.01094 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.005351 on 54 degrees of freedom
## Multiple R-squared: 0.3909, Adjusted R-squared: 0.3007
## F-statistic: 4.332 on 8 and 54 DF, p-value: 0.0004419
testing_ardl7 = ardlDlm(x = scc_test_big$dleaving, y = scc_test_big$d2log_cases, p = 7, q = 7, remove = list(p = c(0,1,2,3,4,5,6), q=c()))
summary(testing_ardl7)
##
## Time series regression with "ts" data:
## Start = 8, End = 70
##
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.060278 -0.004423 0.001805 0.007273 0.048478
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0041398 0.0031885 -1.298 0.19969
## X.7 0.0020389 0.0009014 2.262 0.02774 *
## Y.1 -0.3416279 0.1056659 -3.233 0.00209 **
## Y.2 -0.1038290 0.1125298 -0.923 0.36028
## Y.3 -0.3320478 0.1053011 -3.153 0.00264 **
## Y.4 0.0772674 0.1140214 0.678 0.50088
## Y.5 -0.0936313 0.0996331 -0.940 0.35152
## Y.6 -0.1625195 0.0999608 -1.626 0.10981
## Y.7 0.2141471 0.0787674 2.719 0.00879 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0215 on 54 degrees of freedom
## Multiple R-squared: 0.7712, Adjusted R-squared: 0.7373
## F-statistic: 22.76 on 8 and 54 DF, p-value: 9.715e-15
# all
testing_ardl = ardlDlm(x = scc_test_big$dleaving_mov, y = scc_test_big$d2log_cases_mov, p = 15, q = 15)
summary(testing_ardl)
##
## Time series regression with "ts" data:
## Start = 16, End = 70
##
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.986e-03 -9.923e-04 7.876e-05 1.039e-03 3.087e-03
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0026772 0.0007511 -3.564 0.00165 **
## X.t 0.0016962 0.0011297 1.501 0.14686
## X.1 0.0006484 0.0013834 0.469 0.64368
## X.2 0.0013695 0.0012367 1.107 0.27957
## X.3 -0.0014865 0.0013466 -1.104 0.28105
## X.4 0.0012411 0.0012553 0.989 0.33309
## X.5 0.0002892 0.0011818 0.245 0.80884
## X.6 -0.0002395 0.0011725 -0.204 0.83994
## X.7 0.0008168 0.0012522 0.652 0.52070
## X.8 -0.0004827 0.0013803 -0.350 0.72976
## X.9 0.0004167 0.0011616 0.359 0.72304
## X.10 0.0033635 0.0011412 2.947 0.00723 **
## X.11 -0.0002858 0.0012217 -0.234 0.81713
## X.12 0.0004645 0.0012905 0.360 0.72217
## X.13 -0.0031023 0.0011311 -2.743 0.01159 *
## X.14 0.0035882 0.0011178 3.210 0.00388 **
## X.15 0.0004624 0.0012469 0.371 0.71413
## Y.1 0.0356192 0.1567221 0.227 0.82222
## Y.2 0.0091882 0.1413520 0.065 0.94873
## Y.3 -0.0156067 0.1392002 -0.112 0.91170
## Y.4 -0.0622273 0.1420075 -0.438 0.66533
## Y.5 -0.1140753 0.1147079 -0.994 0.33033
## Y.6 -0.0138435 0.1120046 -0.124 0.90271
## Y.7 -0.3005001 0.1192491 -2.520 0.01913 *
## Y.8 -0.0345115 0.1013978 -0.340 0.73668
## Y.9 -0.0132726 0.0843804 -0.157 0.87639
## Y.10 -0.1900936 0.0851933 -2.231 0.03569 *
## Y.11 -0.0330317 0.0700561 -0.472 0.64172
## Y.12 0.0175546 0.0664028 0.264 0.79385
## Y.13 0.0139956 0.0630154 0.222 0.82620
## Y.14 -0.0201684 0.0639287 -0.315 0.75524
## Y.15 -0.1572236 0.0655611 -2.398 0.02499 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.002179 on 23 degrees of freedom
## Multiple R-squared: 0.9325, Adjusted R-squared: 0.8414
## F-statistic: 10.24 on 31 and 23 DF, p-value: 1.206e-07
testing_ardl = ardlDlm(x = scc_test_big$dleaving, y = scc_test_big$d2log_cases, p = 15, q = 15)
summary(testing_ardl)
##
## Time series regression with "ts" data:
## Start = 16, End = 70
##
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0153804 -0.0046369 0.0002912 0.0038447 0.0200498
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0060194 0.0028369 -2.122 0.044841 *
## X.t 0.0021201 0.0008117 2.612 0.015588 *
## X.1 0.0019892 0.0008597 2.314 0.029966 *
## X.2 0.0035492 0.0009228 3.846 0.000824 ***
## X.3 0.0007635 0.0011123 0.686 0.499316
## X.4 0.0008054 0.0010050 0.801 0.431115
## X.5 0.0006726 0.0009732 0.691 0.496397
## X.6 0.0001779 0.0008291 0.215 0.832030
## X.7 0.0015223 0.0008115 1.876 0.073423 .
## X.8 0.0006827 0.0008229 0.830 0.415259
## X.9 0.0008236 0.0008682 0.949 0.352654
## X.10 0.0021748 0.0008585 2.533 0.018570 *
## X.11 0.0017887 0.0008939 2.001 0.057338 .
## X.12 0.0024315 0.0008839 2.751 0.011385 *
## X.13 0.0001263 0.0009885 0.128 0.899451
## X.14 0.0020863 0.0009276 2.249 0.034374 *
## X.15 0.0006511 0.0009032 0.721 0.478274
## Y.1 -0.7371592 0.2006844 -3.673 0.001262 **
## Y.2 -0.5358158 0.1785719 -3.001 0.006382 **
## Y.3 -0.3932115 0.1892215 -2.078 0.049052 *
## Y.4 -0.2606945 0.1764794 -1.477 0.153186
## Y.5 -0.2385887 0.1541623 -1.548 0.135358
## Y.6 -0.0653172 0.1382773 -0.472 0.641119
## Y.7 -0.1650019 0.1254548 -1.315 0.201395
## Y.8 -0.2115432 0.1224248 -1.728 0.097404 .
## Y.9 -0.2025611 0.1114602 -1.817 0.082218 .
## Y.10 -0.3257234 0.1151342 -2.829 0.009512 **
## Y.11 -0.2915754 0.1271660 -2.293 0.031332 *
## Y.12 -0.1639630 0.1302824 -1.259 0.220829
## Y.13 -0.0242918 0.1241679 -0.196 0.846613
## Y.14 0.0824732 0.1067281 0.773 0.447545
## Y.15 0.0886663 0.0772121 1.148 0.262629
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01069 on 23 degrees of freedom
## Multiple R-squared: 0.8363, Adjusted R-squared: 0.6156
## F-statistic: 3.79 on 31 and 23 DF, p-value: 0.0007732
# try with not log of cases
testing_ardl_not_log = ardlDlm(x = scc_test_big$dleaving, y = scc_test_big$d2cases, p = 15, q = 15)
summary(testing_ardl_not_log)
##
## Time series regression with "ts" data:
## Start = 16, End = 70
##
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.8077 -4.6345 0.3506 5.5457 14.4709
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.03807 1.87360 0.020 0.983965
## X.t 1.06282 0.89591 1.186 0.247616
## X.1 1.72185 0.84293 2.043 0.052708 .
## X.2 2.08779 0.76098 2.744 0.011574 *
## X.3 -0.29649 0.91080 -0.326 0.747724
## X.4 -1.03824 1.00674 -1.031 0.313130
## X.5 -1.98815 0.95116 -2.090 0.047845 *
## X.6 -0.29275 0.84875 -0.345 0.733294
## X.7 0.77840 0.85529 0.910 0.372208
## X.8 0.58843 0.81083 0.726 0.475340
## X.9 -0.37824 0.76945 -0.492 0.627685
## X.10 0.31816 0.71670 0.444 0.661249
## X.11 0.27307 0.68208 0.400 0.692594
## X.12 0.82349 0.66682 1.235 0.229316
## X.13 -0.74080 0.68498 -1.081 0.290689
## X.14 2.53881 0.76084 3.337 0.002864 **
## X.15 0.21156 0.82757 0.256 0.800505
## Y.1 -1.47509 0.18745 -7.869 5.68e-08 ***
## Y.2 -1.70247 0.28862 -5.899 5.18e-06 ***
## Y.3 -1.83794 0.36017 -5.103 3.62e-05 ***
## Y.4 -1.94224 0.42495 -4.571 0.000136 ***
## Y.5 -1.73150 0.46087 -3.757 0.001026 **
## Y.6 -1.62403 0.46350 -3.504 0.001910 **
## Y.7 -1.51701 0.46735 -3.246 0.003564 **
## Y.8 -1.47255 0.46274 -3.182 0.004151 **
## Y.9 -1.33424 0.45187 -2.953 0.007139 **
## Y.10 -1.25388 0.42155 -2.974 0.006785 **
## Y.11 -1.04299 0.38103 -2.737 0.011740 *
## Y.12 -0.85170 0.34327 -2.481 0.020842 *
## Y.13 -0.58238 0.29787 -1.955 0.062827 .
## Y.14 -0.21879 0.22868 -0.957 0.348633
## Y.15 -0.10857 0.14319 -0.758 0.455987
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.23 on 23 degrees of freedom
## Multiple R-squared: 0.9325, Adjusted R-squared: 0.8414
## F-statistic: 10.24 on 31 and 23 DF, p-value: 1.207e-07
# pre april 15th
scc_test_pre_april15 <- scc_test_big %>%
filter(date <= "2020-04-15")
# all with this dataset
testing_ardl_pre = ardlDlm(x = scc_test_pre_april15$dleaving_mov, y = scc_test_pre_april15$d2log_cases_mov, p = 4, q = 4, remove = list(p = c(0,1,2,3)))
summary(testing_ardl_pre)
##
## Time series regression with "ts" data:
## Start = 5, End = 46
##
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0193970 -0.0027041 -0.0002263 0.0023051 0.0263998
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.002201 0.001787 -1.231 0.22617
## X.4 -0.001325 0.001954 -0.678 0.50215
## Y.1 -0.053757 0.168624 -0.319 0.75172
## Y.2 0.248812 0.176541 1.409 0.16731
## Y.3 0.439239 0.143173 3.068 0.00408 **
## Y.4 0.069857 0.155962 0.448 0.65690
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.009699 on 36 degrees of freedom
## Multiple R-squared: 0.2384, Adjusted R-squared: 0.1327
## F-statistic: 2.254 on 5 and 36 DF, p-value: 0.06978
testing_ardl_pre = ardlDlm(x = scc_test_pre_april15$dleaving, y = scc_test_pre_april15$d2log_cases, p = 4, q = 4, remove = list(p = c(0,1,2,3)))
summary(testing_ardl_pre)
##
## Time series regression with "ts" data:
## Start = 5, End = 46
##
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.101504 -0.022735 -0.004123 0.018615 0.098219
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.009738 0.007404 -1.315 0.196778
## X.4 0.002415 0.001969 1.227 0.227944
## Y.1 -0.697380 0.162307 -4.297 0.000126 ***
## Y.2 -0.597481 0.179501 -3.329 0.002022 **
## Y.3 -0.469015 0.187412 -2.503 0.017008 *
## Y.4 -0.213756 0.138632 -1.542 0.131844
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04614 on 36 degrees of freedom
## Multiple R-squared: 0.3643, Adjusted R-squared: 0.276
## F-statistic: 4.127 on 5 and 36 DF, p-value: 0.004587
# ardlBoundOrders
ardlBoundOrders(data = scc_test_big %>% dplyr::select(dleaving, d2log_cases), d2log_cases~dleaving)
## $p
## dleaving
## 1 14
##
## $q
## [1] 15
##
## $Stat.table
## q = 1 q = 2 q = 3 q = 4 q = 5 q = 6 q = 7
## p = 1 -250.9489 -267.5710 -263.3483 -275.1213 -293.5316 -305.3056 -297.1058
## p = 2 -256.9344 -267.4827 -265.8688 -274.3035 -294.2032 -307.1245 -299.1835
## p = 3 -263.4581 -263.4581 -268.7743 -273.5211 -292.4649 -305.1729 -297.3040
## p = 4 -256.9197 -272.2833 -272.2833 -275.1706 -290.9673 -303.4811 -295.7392
## p = 5 -258.4962 -280.4359 -279.0057 -279.0057 -288.9916 -301.9384 -294.4630
## p = 6 -280.4465 -298.7995 -296.9824 -297.1321 -297.1321 -302.8457 -295.1018
## p = 7 -271.7965 -299.8338 -298.2299 -298.3701 -297.4849 -297.4849 -296.7258
## p = 8 -282.1093 -292.7891 -290.7999 -288.9308 -295.7467 -295.3447 -295.3447
## p = 9 -279.1848 -291.0854 -289.6336 -287.9535 -291.7501 -291.7010 -295.5427
## p = 10 -294.6798 -298.3682 -296.6957 -297.1293 -297.0726 -295.1486 -294.7634
## p = 11 -271.2616 -280.4542 -278.5015 -276.5136 -281.9553 -280.0605 -280.1146
## p = 12 -289.7541 -293.7644 -294.2128 -296.9428 -295.2164 -293.2390 -291.4248
## p = 13 -310.9712 -314.3384 -318.1372 -320.8433 -319.1747 -322.3125 -320.3157
## p = 14 -305.1308 -310.9983 -311.7019 -320.4972 -319.1592 -318.0871 -317.2385
## p = 15 -324.2096 -323.0515 -328.8318 -329.7601 -328.1312 -327.9968 -326.5103
## q = 8 q = 9 q = 10 q = 11 q = 12 q = 13 q = 14
## p = 1 -304.7931 -297.7748 -293.8927 -290.5643 -295.6733 -293.3990 -298.1992
## p = 2 -306.1365 -299.4384 -293.9848 -290.1470 -294.4886 -292.6269 -296.7515
## p = 3 -305.0713 -298.0177 -293.0737 -288.9128 -294.2433 -292.4169 -295.4953
## p = 4 -303.5975 -297.1456 -291.9681 -287.3530 -292.3349 -290.5310 -294.4090
## p = 5 -304.7066 -299.1953 -292.0180 -288.2497 -291.5839 -289.3112 -294.7581
## p = 6 -303.4262 -298.3603 -291.5864 -287.5411 -291.7405 -290.0525 -297.9255
## p = 7 -302.1107 -297.2501 -289.9498 -285.5484 -290.5526 -291.0481 -298.6609
## p = 8 -303.6832 -297.1590 -290.5779 -284.6231 -288.6746 -289.3063 -296.7745
## p = 9 -295.5427 -297.4785 -291.6063 -286.6418 -288.2629 -292.4720 -297.9880
## p = 10 -293.5849 -293.5849 -293.6999 -292.4071 -291.5141 -293.3995 -301.3351
## p = 11 -283.8958 -285.8242 -285.8242 -296.8939 -301.7199 -303.3599 -318.2587
## p = 12 -290.6325 -289.5466 -300.0136 -300.0136 -299.7320 -304.6793 -318.1959
## p = 13 -318.3452 -320.1624 -321.2812 -329.1462 -329.1462 -327.7065 -325.9072
## p = 14 -315.3132 -317.4421 -318.5133 -324.4961 -324.0695 -324.0695 -325.1359
## p = 15 -325.0579 -330.1588 -329.0686 -333.6106 -340.9263 -341.0082 -341.0082
## q = 15
## p = 1 -317.9986
## p = 2 -316.0589
## p = 3 -314.5318
## p = 4 -312.8589
## p = 5 -314.6978
## p = 6 -315.5104
## p = 7 -320.9857
## p = 8 -318.9927
## p = 9 -317.3532
## p = 10 -320.8633
## p = 11 -330.8053
## p = 12 -329.2307
## p = 13 -342.8803
## p = 14 -343.2450
## p = 15 -342.2451
##
## $min.Stat
## [1] -343.245
ARDL test might not require stationarity?
# 4th order lag on x only, original variables
testing_ardl4_2 = ardlDlm(x = scc_test_big$leaving_home_dif, y = scc_test_big$cases, p = 4, q = 4, remove = list(p = c(0,1,2,3), q=c(1,2,3,4)))
summary(testing_ardl4_2)
##
## Time series regression with "ts" data:
## Start = 5, End = 70
##
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.619e-14 -1.373e-15 -3.200e-16 5.980e-16 5.774e-14
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.653e-15 2.920e-15 2.279e+00 0.0261 *
## X.4 -2.633e-15 2.065e-16 -1.275e+01 < 2e-16 ***
## Y.1 4.761e-16 6.643e-17 7.167e+00 1.11e-09 ***
## Y.0 1.000e+00 6.854e-17 1.459e+16 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.974e-15 on 62 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 2.166e+35 on 3 and 62 DF, p-value: < 2.2e-16
# ardlBoundOrders
ardlBoundOrders(data = scc_test_big %>% dplyr::select(leaving_home_dif, cases), cases~leaving_home_dif)
## $p
## leaving_home_dif
## 1 15
##
## $q
## [1] 5
##
## $Stat.table
## q = 1 q = 2 q = 3 q = 4 q = 5 q = 6 q = 7 q = 8
## p = 1 551.4315 541.7830 535.8696 530.0526 524.9652 519.4883 505.4713 495.0153
## p = 2 541.2776 542.2030 534.6996 528.1887 523.1593 517.9178 505.7249 495.1092
## p = 3 534.7343 534.7343 536.6990 530.0122 524.8594 519.7008 507.5232 496.6103
## p = 4 528.1732 529.7440 529.7440 530.8711 524.8219 519.0623 507.7080 495.9187
## p = 5 521.6816 522.9820 524.8233 524.8233 526.8219 520.8457 508.9555 495.2984
## p = 6 512.8850 514.6301 516.4811 518.0836 518.0836 519.5213 509.8541 496.7291
## p = 7 506.1375 508.0238 509.6129 511.0533 512.9293 512.9293 511.4407 497.9630
## p = 8 498.6552 500.3424 501.2760 501.6296 502.9350 503.2506 503.2506 495.5341
## p = 9 485.4073 486.0239 486.9155 487.9907 489.8081 491.6142 482.9317 482.9317
## p = 10 476.5659 476.4151 477.7322 477.6431 479.4408 480.9793 476.9776 476.8845
## p = 11 462.1229 462.0401 463.9358 464.6950 466.6923 467.5709 463.7286 463.3197
## p = 12 452.4190 449.4032 451.3701 452.5185 454.4904 456.3255 455.2871 454.5650
## p = 13 446.0101 443.5004 444.7262 441.0506 443.0246 444.7632 446.1498 445.3551
## p = 14 414.3307 413.6259 414.6902 415.8575 414.6262 414.8230 416.5557 417.7452
## p = 15 412.1510 406.2678 406.6905 405.6631 404.3863 405.3680 406.8105 407.5832
## q = 9 q = 10 q = 11 q = 12 q = 13 q = 14 q = 15
## p = 1 487.5733 482.6726 474.7761 469.6282 463.6572 455.9885 446.9640
## p = 2 487.4514 481.4059 472.0640 466.7745 460.9846 452.7444 446.1624
## p = 3 489.1777 483.1421 472.1753 466.7051 460.7315 452.3941 446.6445
## p = 4 489.3430 483.7567 473.7310 468.1435 462.0371 452.7884 447.6725
## p = 5 488.0198 482.1261 474.4132 469.3686 463.1432 452.8999 447.4158
## p = 6 489.9574 484.1184 475.9503 470.3785 463.1308 453.2783 448.0971
## p = 7 490.3463 483.2657 474.3034 469.7507 463.7878 454.7004 449.0435
## p = 8 484.6181 475.1708 462.4046 457.7744 452.6782 444.1399 438.7220
## p = 9 484.9306 476.5633 464.3612 459.7471 454.5246 446.0604 440.2099
## p = 10 476.8845 477.5860 466.1511 461.6051 456.0736 446.5651 440.2906
## p = 11 464.3016 464.3016 466.2455 461.6123 456.7817 448.4575 441.8958
## p = 12 456.2021 457.6613 457.6613 459.4294 454.9623 448.3276 439.1265
## p = 13 447.1687 448.1837 450.1686 450.1686 452.0140 443.3848 436.2394
## p = 14 419.6475 420.7660 422.2259 423.2039 423.2039 425.1985 410.9562
## p = 15 409.1554 409.5049 411.4599 411.3893 413.3596 413.3596 409.0664
##
## $min.Stat
## [1] 404.3863
ardlBoundOrders(data = scc_test_big %>% dplyr::select(leaving_home_dif, dcases), dcases~leaving_home_dif)
## $p
## leaving_home_dif
## 1 15
##
## $q
## [1] 4
##
## $Stat.table
## q = 1 q = 2 q = 3 q = 4 q = 5 q = 6 q = 7 q = 8
## p = 1 563.3203 556.5468 550.6648 543.6184 536.7660 525.5309 512.5341 507.2811
## p = 2 556.9574 558.4840 552.6581 545.6140 538.7639 527.3610 514.4775 509.2262
## p = 3 549.4072 549.4072 551.3550 545.3277 538.0811 526.3433 515.1568 509.6884
## p = 4 543.8399 545.2201 545.2201 546.6075 539.6319 527.2861 516.5984 511.1939
## p = 5 535.6211 536.9059 538.8974 538.8974 539.7245 527.6348 517.8986 512.6875
## p = 6 521.7455 523.4562 525.2409 526.6596 526.6596 523.3985 513.2089 508.4337
## p = 7 511.8901 513.2635 515.1744 516.7433 518.4597 518.4597 514.8977 510.1451
## p = 8 515.9819 516.6476 517.4174 516.9702 516.4896 510.0367 510.0367 511.9807
## p = 9 495.7937 497.7929 498.8826 500.0176 501.4963 495.5237 493.2470 493.2470
## p = 10 486.0202 488.0184 489.8807 490.1027 491.1960 484.7809 485.1142 485.6264
## p = 11 476.8587 477.9075 479.9013 480.1035 481.7970 476.3355 475.3342 475.8643
## p = 12 466.2053 466.9407 467.8379 469.6922 471.1108 467.6380 467.1936 468.2445
## p = 13 448.7376 450.7376 452.4398 450.8206 452.6976 451.3775 452.9695 452.4121
## p = 14 430.0209 431.8791 433.0050 423.6393 425.2960 426.3064 428.3032 428.3932
## p = 15 423.2949 425.1817 427.0081 421.2607 423.2438 421.3601 423.3568 421.9097
## q = 9 q = 10 q = 11 q = 12 q = 13 q = 14 q = 15
## p = 1 500.9525 495.3731 488.5888 483.6054 478.3645 470.1307 462.2284
## p = 2 502.8943 497.3640 490.5714 485.6030 480.3645 472.0792 464.2245
## p = 3 503.0829 497.3656 491.5017 486.5515 481.0187 473.0953 464.8591
## p = 4 503.9822 496.9930 490.0851 484.9604 479.5173 472.2849 463.1397
## p = 5 505.6865 498.1247 489.4702 484.4219 478.6026 470.9423 463.6693
## p = 6 502.7776 495.0681 483.5585 476.0396 469.9647 461.4657 454.9348
## p = 7 504.6644 497.0398 485.5517 477.8507 471.0057 463.2649 456.8309
## p = 8 506.3176 498.7452 486.2678 478.4237 472.3820 461.7054 454.7676
## p = 9 495.1670 488.6104 476.0440 469.5490 461.5124 454.3281 446.5764
## p = 10 485.6264 485.9962 474.3965 466.8136 457.1185 450.2643 442.2470
## p = 11 477.8021 477.8021 476.3958 468.8136 458.5311 451.7515 444.2388
## p = 12 470.2362 471.3915 471.3915 470.2410 460.2490 453.2161 444.4992
## p = 13 454.4119 454.3142 455.4411 455.4411 451.4153 444.5477 436.7992
## p = 14 430.3769 430.0190 431.9277 432.4159 432.4159 432.6255 422.5280
## p = 15 423.8915 423.3048 425.1451 423.5634 423.5652 423.5652 424.4842
##
## $min.Stat
## [1] 421.2607
ardlBoundOrders(data = scc_test_big %>% dplyr::select(leaving_home_dif, d2cases), d2cases~leaving_home_dif)
## $p
## leaving_home_dif
## 1 14
##
## $q
## [1] 15
##
## $Stat.table
## q = 1 q = 2 q = 3 q = 4 q = 5 q = 6 q = 7 q = 8
## p = 1 596.8029 585.1344 574.8355 567.6410 534.8172 527.7619 520.8212 515.8249
## p = 2 590.1137 587.0190 576.8176 569.6169 536.7224 529.5562 522.6747 517.6850
## p = 3 576.8215 576.8215 575.6664 568.6428 537.9695 530.9054 524.3234 519.3164
## p = 4 575.9275 569.5017 569.5017 570.2982 539.7050 532.7799 526.2488 521.2463
## p = 5 570.1679 566.1116 564.1554 564.1554 541.6558 534.6894 528.2101 523.2170
## p = 6 558.7235 554.8029 554.6769 534.3003 534.3003 534.5655 527.6211 522.5500
## p = 7 546.1128 542.9943 540.9415 539.9070 528.8434 528.8434 529.0404 524.0259
## p = 8 548.2126 543.6979 540.9000 541.3655 522.3490 522.3181 522.3181 524.2777
## p = 9 535.6439 530.0268 526.6979 528.2749 510.9973 512.6614 513.7313 513.7313
## p = 10 529.4588 524.6606 521.2560 522.7175 506.5407 508.4563 508.6884 510.4387
## p = 11 517.4364 513.5564 508.7137 510.1027 493.8151 495.5684 496.0371 498.0363
## p = 12 506.4488 502.5069 499.8967 501.0570 484.9590 486.6526 485.8424 487.8103
## p = 13 481.8881 480.3362 479.9163 478.9870 471.8291 473.7712 473.8086 475.6352
## p = 14 463.0507 463.3476 454.9262 450.7764 448.0149 448.6843 445.7275 447.4288
## p = 15 445.4927 446.3352 438.4082 440.4080 437.1953 438.5188 436.6688 438.3851
## q = 9 q = 10 q = 11 q = 12 q = 13 q = 14 q = 15
## p = 1 509.7835 504.5044 498.9399 486.8007 479.8890 473.8559 467.4932
## p = 2 511.6966 506.4000 500.7271 488.6346 481.2843 475.3213 469.2805
## p = 3 513.1026 507.7074 502.3867 490.1641 482.5489 476.9465 470.6609
## p = 4 514.8333 509.3236 503.7799 491.3077 483.1360 476.8671 471.0562
## p = 5 516.8324 511.3215 505.7229 492.9721 484.5546 477.8199 472.5858
## p = 6 516.8858 511.6751 505.8550 492.9150 484.2152 476.3484 471.2303
## p = 7 518.0338 512.7477 506.2194 492.8402 484.1382 475.5451 470.3127
## p = 8 518.0536 512.3625 504.3034 488.8966 478.0287 464.8852 460.1567
## p = 9 515.3571 509.4286 502.4835 487.2160 475.5037 463.7814 458.4289
## p = 10 510.4387 511.3092 504.2173 488.8262 477.4583 465.1542 459.7079
## p = 11 499.4619 499.4619 500.6853 485.8143 475.4783 464.9744 458.4922
## p = 12 489.5813 490.9129 490.9129 480.0151 469.4743 460.5481 452.5700
## p = 13 477.3555 478.8883 470.2526 470.2526 469.8552 460.5796 452.6569
## p = 14 446.0042 446.6970 444.8195 439.8721 439.8721 439.9391 430.7741
## p = 15 436.0841 437.8856 438.1497 431.7956 432.0592 432.0592 432.7404
##
## $min.Stat
## [1] 430.7741
ardlBoundOrders(data = scc_test_big %>% dplyr::select(leaving_home_dif, d2log_cases), d2log_cases~leaving_home_dif)
## $p
## leaving_home_dif
## 1 15
##
## $q
## [1] 15
##
## $Stat.table
## q = 1 q = 2 q = 3 q = 4 q = 5 q = 6 q = 7
## p = 1 -247.1158 -262.6566 -259.7280 -271.1864 -289.5665 -306.0241 -297.9282
## p = 2 -260.1151 -270.6964 -271.7295 -276.1488 -293.3896 -305.4478 -297.7216
## p = 3 -272.1176 -272.1176 -275.2289 -276.2051 -295.3421 -308.4182 -300.1463
## p = 4 -257.3371 -276.7596 -276.7596 -275.8190 -293.9966 -306.7898 -298.4602
## p = 5 -257.5284 -277.3622 -275.5188 -275.5188 -292.2041 -304.8192 -296.4861
## p = 6 -273.3595 -291.6139 -290.3658 -291.6970 -291.6970 -303.0178 -294.7003
## p = 7 -272.3023 -300.7645 -299.2209 -297.8259 -297.7265 -297.7265 -296.0480
## p = 8 -282.0520 -292.6506 -290.6616 -288.8317 -296.8763 -295.5956 -295.5956
## p = 9 -278.5346 -293.5507 -292.1824 -290.5130 -294.2584 -292.7547 -295.4100
## p = 10 -292.4685 -297.2686 -295.9552 -295.7238 -295.5189 -293.5769 -292.4145
## p = 11 -271.5680 -281.0769 -279.1167 -277.1916 -282.6667 -280.7441 -280.4900
## p = 12 -297.7511 -299.3264 -302.0489 -303.4335 -301.4789 -299.6070 -297.8830
## p = 13 -294.5457 -294.4147 -296.4900 -301.1529 -299.8009 -299.7930 -299.7014
## p = 14 -306.9080 -314.9017 -318.8187 -328.3627 -326.4117 -325.5119 -323.5907
## p = 15 -322.5124 -320.8542 -327.1413 -330.1459 -329.1888 -328.1993 -326.7715
## q = 8 q = 9 q = 10 q = 11 q = 12 q = 13 q = 14
## p = 1 -303.1134 -296.7856 -293.6299 -289.3689 -299.7575 -300.6864 -295.6658
## p = 2 -302.9945 -295.8307 -292.4216 -288.5985 -301.3512 -307.0429 -306.3532
## p = 3 -304.6415 -297.5186 -292.3243 -288.1629 -299.7060 -306.6731 -307.4808
## p = 4 -303.5515 -296.1163 -291.4783 -286.9199 -299.7326 -306.3824 -308.1563
## p = 5 -302.3817 -295.2922 -290.2753 -285.3606 -297.7328 -304.5140 -306.6937
## p = 6 -302.8127 -297.2628 -290.6908 -286.4765 -298.3092 -305.8850 -308.5185
## p = 7 -301.7097 -296.3767 -289.7787 -285.6182 -297.0538 -304.7369 -307.9607
## p = 8 -300.4574 -295.3230 -288.0325 -283.6450 -295.5137 -305.0929 -309.1502
## p = 9 -295.4100 -295.2668 -288.7073 -282.8341 -294.0422 -303.4526 -307.4371
## p = 10 -291.3980 -291.3980 -289.6450 -284.6432 -292.1702 -303.2680 -306.9536
## p = 11 -283.4290 -285.4861 -285.4861 -290.4073 -294.3568 -302.0437 -309.2631
## p = 12 -297.7796 -295.7838 -302.8100 -302.8100 -305.2507 -312.4671 -322.4358
## p = 13 -298.3913 -296.9694 -298.4937 -302.1586 -302.1586 -314.0261 -322.6509
## p = 14 -322.4463 -325.2010 -324.5372 -329.7972 -328.1825 -328.1825 -331.6931
## p = 15 -325.3757 -330.3256 -329.2409 -333.4078 -340.3607 -341.1173 -341.1173
## q = 15
## p = 1 -311.9700
## p = 2 -317.2822
## p = 3 -315.7170
## p = 4 -313.7580
## p = 5 -311.9534
## p = 6 -314.2362
## p = 7 -314.2484
## p = 8 -319.3696
## p = 9 -317.3698
## p = 10 -315.6900
## p = 11 -319.0116
## p = 12 -328.8293
## p = 13 -327.2857
## p = 14 -340.9265
## p = 15 -341.4775
##
## $min.Stat
## [1] -341.4775
ardlBoundOrders(data = scc_test_big %>% dplyr::select(leaving_home_dif, cases_growth_daily), cases_growth_daily~leaving_home_dif)
## $p
## leaving_home_dif
## 1 15
##
## $q
## [1] 15
##
## $Stat.table
## q = 1 q = 2 q = 3 q = 4 q = 5 q = 6 q = 7 q = 8
## p = 1 393.1757 375.3039 357.0862 340.7902 334.0467 282.5013 275.6036 269.8120
## p = 2 379.1824 376.6348 352.4453 339.6024 331.6322 284.2639 277.5686 271.7322
## p = 3 369.7678 369.7678 353.5871 340.5742 332.7598 285.2172 277.9381 272.7965
## p = 4 341.3875 343.1355 343.1355 340.2579 334.4932 286.7561 279.8589 274.3915
## p = 5 347.9692 347.9300 334.7346 334.7346 329.4404 288.4827 281.6367 276.3697
## p = 6 331.4535 332.4995 313.9106 312.9811 312.9811 290.0421 282.5161 276.0033
## p = 7 296.1998 295.4014 282.5097 284.3730 284.3179 284.3179 284.0541 277.7899
## p = 8 289.1624 290.9319 285.0433 286.7694 288.5234 276.8694 276.8694 278.8337
## p = 9 289.7759 291.0523 276.4356 278.4331 279.7860 271.4070 273.2216 273.2216
## p = 10 276.5552 278.0002 266.3790 267.5950 267.8485 257.1157 259.0252 260.6031
## p = 11 262.4497 263.5248 261.6312 262.3240 264.3200 261.7908 263.7794 264.2309
## p = 12 221.5408 221.4510 222.8320 223.6527 221.4774 223.4773 223.8759 225.7048
## p = 13 209.1545 206.8361 207.7154 208.5987 202.8317 204.6966 201.2008 198.7951
## p = 14 198.2758 198.5349 189.7918 185.3501 175.6170 177.6005 179.4515 180.1959
## p = 15 191.6988 185.1561 186.7437 180.1345 175.2661 175.9340 176.3101 177.5464
## q = 9 q = 10 q = 11 q = 12 q = 13 q = 14 q = 15
## p = 1 264.4046 256.5817 251.1818 218.0678 203.8873 198.1708 182.5250
## p = 2 266.4033 258.4472 253.0561 219.7598 204.2247 197.6611 180.6552
## p = 3 266.7027 259.5050 254.8955 221.6009 206.1549 198.7904 181.0040
## p = 4 268.4987 261.2078 256.5957 223.2449 208.1431 200.2166 182.9999
## p = 5 270.1325 262.8564 258.5765 224.8884 209.9436 201.5634 184.0590
## p = 6 269.1356 260.2955 256.8728 220.7890 199.4582 192.4838 175.7583
## p = 7 270.9138 261.6299 258.8486 222.7881 200.2446 190.4380 175.9913
## p = 8 272.7483 263.6298 260.8486 224.7564 202.1607 192.4373 176.3502
## p = 9 267.5247 259.9905 258.4628 225.1660 200.5352 189.2311 177.3842
## p = 10 260.6031 261.2215 258.7955 226.5752 201.8393 190.2866 179.1897
## p = 11 260.3329 260.3329 259.4055 227.6676 203.3513 191.5684 180.2022
## p = 12 224.9154 226.8497 226.8497 225.0855 194.6972 188.4177 171.0937
## p = 13 197.6100 199.4856 199.3437 199.3437 196.3416 188.1092 172.8223
## p = 14 170.8064 171.9388 173.2853 168.8885 168.8885 170.6939 160.8284
## p = 15 179.4827 177.8003 175.8217 170.5812 159.8029 159.8029 161.4884
##
## $min.Stat
## [1] 159.8029
testing_ardl4_growth = ardlDlm(x = scc_test_big$leaving_home_dif, y = scc_test_big$cases_growth_daily, p = 5, q = 5, remove = list(p = c(0,1,2,3)))
summary(testing_ardl4_growth)
##
## Time series regression with "ts" data:
## Start = 6, End = 70
##
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.7931 -1.1739 -0.0644 1.1413 12.4465
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.77637 4.04403 0.687 0.495158
## X.4 0.28624 0.15813 1.810 0.075547 .
## X.5 -0.19359 0.16047 -1.206 0.232641
## Y.1 0.23420 0.11525 2.032 0.046815 *
## Y.2 0.19142 0.11495 1.665 0.101355
## Y.3 -0.09892 0.11362 -0.871 0.387578
## Y.4 0.37583 0.09750 3.855 0.000297 ***
## Y.5 0.08068 0.09769 0.826 0.412309
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.755 on 57 degrees of freedom
## Multiple R-squared: 0.8358, Adjusted R-squared: 0.8156
## F-statistic: 41.43 on 7 and 57 DF, p-value: < 2.2e-16
scc_nolag <- scc_test_big %>%
dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving_mov) %>%
cbind(`Lag` = "0")
scc_3lag <- scc_test_big %>%
dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving3_mov) %>%
cbind(`Lag` = "3")
colnames(scc_3lag)[ncol(scc_3lag)-1] = "leaving_mov"
scc_4lag <- scc_test_big %>%
dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving4_mov) %>%
cbind(`Lag` = "4")
colnames(scc_4lag)[ncol(scc_4lag)-1] = "leaving_mov"
scc_5lag <- scc_test_big %>%
dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving5_mov) %>%
cbind(`Lag` = "5")
colnames(scc_5lag)[ncol(scc_5lag)-1] = "leaving_mov"
scc_6lag <- scc_test_big %>%
dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving6_mov) %>%
cbind(`Lag` = "6")
colnames(scc_6lag)[ncol(scc_6lag)-1] = "leaving_mov"
scc_7lag <- scc_test_big %>%
dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving7_mov) %>%
cbind(`Lag` = "7")
colnames(scc_7lag)[ncol(scc_7lag)-1] = "leaving_mov"
scc_8lag <- scc_test_big %>%
dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving8_mov) %>%
cbind(`Lag` = "8")
colnames(scc_8lag)[ncol(scc_8lag)-1] = "leaving_mov"
scc_9lag <- scc_test_big %>%
dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving9_mov) %>%
cbind(`Lag` = "9")
colnames(scc_9lag)[ncol(scc_9lag)-1] = "leaving_mov"
scc_10lag <- scc_test_big %>%
dplyr::select(date, dcases, cases_growth_daily_mov, cases, log_cases, leaving10_mov) %>%
cbind(`Lag` = "10")
colnames(scc_10lag)[ncol(scc_10lag)-1] = "leaving_mov"
scc_scatter <- rbind(scc_nolag, scc_3lag, scc_4lag, scc_5lag, scc_6lag, scc_7lag, scc_8lag, scc_9lag, scc_10lag)
fig_cases_growth <-
plot_ly (scc_scatter) %>%
add_trace(
x = ~leaving_mov,
y = ~cases_growth_daily_mov,
frame = ~`Lag`,
type = 'scatter',
mode = 'markers',
showlegend = F
) %>%
animation_button(visible = F) %>%
layout(xaxis = list(title = '% difference in leaving home'), yaxis = list(title = 'daily case growth rate'), title = list(title = "SCC"))
fig_cases_growth
fig_cases_abs <-
plot_ly (scc_scatter) %>%
add_trace(
x = ~leaving_mov,
y = ~dcases,
frame = ~`Lag`,
type = 'scatter',
mode = 'markers',
showlegend = F
) %>%
animation_button(visible = F) %>%
layout(xaxis = list(title = '% difference in leaving home'), yaxis = list(title = 'daily new cases'), title = list(title = "SCC"))
fig_cases_abs
For now, just plotting the metric of cumulative percent positive
smc_cases_sd_daily <- bay_sd %>%
filter(origin_census_block_group %in% smc_blockgroups$GEOID) %>%
group_by(date) %>%
summarize(total_at_home = sum(completely_home_device_count), total_devices = sum(device_count)) %>%
mutate(
percent_at_home = total_at_home*100/total_devices,
percent_leaving_home = (100 - percent_at_home),
) %>%
left_join(
smc_cumulative_cases
) %>%
filter(date >= (min(smc_cumulative_cases$date)-10))
# get the baseline percent of people at home
pre_case_growth_at_home_smc <- bay_sd %>%
filter(date < min(smc_cumulative_cases$date)) %>%
filter(origin_census_block_group %in% smc_blockgroups$GEOID) %>%
summarize(total_at_home = sum(completely_home_device_count), total_devices = sum(device_count)) %>%
mutate(percent_at_home = total_at_home*100/total_devices, percent_leaving_home = (100 - percent_at_home))
# include change in percent leaving home
smc_cases_sd_daily <- smc_cases_sd_daily %>%
mutate(leaving_home_dif = percent_leaving_home - pre_case_growth_at_home_smc$percent_leaving_home[1],
log_cases = log(cumulative_cases))
smc_test_small <- smc_cases_sd_daily %>%
mutate(
leaving_mov = movavg(leaving_home_dif, 7, type = "s"),
leaving4 = c(rep(NA,4), smc_cases_sd_daily$leaving_home_dif[1:(nrow(smc_cases_sd_daily)-4)]),
leaving4_mov = movavg(leaving4, 7, type = "s"),
leaving5 = c(rep(NA,5), smc_cases_sd_daily$leaving_home_dif[1:(nrow(smc_cases_sd_daily)-5)]),
leaving5_mov = movavg(leaving5, 7, type = "s"),
leaving6 = c(rep(NA,6), smc_cases_sd_daily$leaving_home_dif[1:(nrow(smc_cases_sd_daily)-6)]),
leaving6_mov = movavg(leaving6, 7, type = "s"),
leaving7 = c(rep(NA,7), smc_cases_sd_daily$leaving_home_dif[1:(nrow(smc_cases_sd_daily)-7)]),
leaving7_mov = movavg(leaving7, 7, type = "s"),
leaving8 = c(rep(NA,8), smc_cases_sd_daily$leaving_home_dif[1:(nrow(smc_cases_sd_daily)-8)]),
leaving8_mov = movavg(leaving8, 7, type = "s"),
leaving9 = c(rep(NA,9), smc_cases_sd_daily$leaving_home_dif[1:(nrow(smc_cases_sd_daily)-9)]),
leaving9_mov = movavg(leaving9, 7, type = "s"),
dcases = c(NA,diff(cumulative_cases)),
leaving10 = c(rep(NA,10), smc_cases_sd_daily$leaving_home_dif[1:(nrow(smc_cases_sd_daily)-10)]),
leaving10_mov = movavg(leaving10, 7, type = "s"),
dcases = c(NA,diff(cumulative_cases)),
past_cases = c(NA, smc_cases_sd_daily$cumulative_cases[1:(nrow(smc_cases_sd_daily)-1)]),
cases_growth_daily = (dcases / past_cases)*100,
cases_growth_daily_mov = movavg(cases_growth_daily, 7, type = "s")) %>%
filter(date >= "2020-03-01")
smc_test_small %>% ggplot(
aes(x = date)) +
geom_line(aes(y = leaving6, color="Leaving home")) +
geom_line(aes(y = cases_growth_daily-30, color = "Cases")) +
scale_y_continuous(sec.axis = sec_axis(~.*1+30, name = "Daily case growth rate (%)")) +
scale_colour_manual(values = c("red", "blue")) +
labs(y = "Change in percent of devices leaving home 6 days ago relative to before", x = "Date", color = "Data", title = "San Mateo County Case Growth and Social Distancing, 6 Day Lag")
# moving average
smc_test_small %>% ggplot(
aes(x = date)) +
geom_line(aes(y = leaving6_mov, color="Leaving home")) +
geom_line(aes(y = cases_growth_daily_mov-30, color = "Cases")) +
scale_y_continuous(sec.axis = sec_axis(~.*1+30, name = "Daily case growth rate (%), 7 day moving average")) +
scale_colour_manual(values = c("red", "blue")) +
labs(y = "Change in percent of devices leaving home 6 days ago relative to before, 7 day moving average", x = "Date", color = "Data", title = "San Mateo County Case Growth and Social Distancing, 6 Day Lag")
# try with percent positive cumulatively
smc_test_small %>% ggplot(
aes(x = date)) +
geom_line(aes(y = leaving6, color="Leaving home")) +
geom_line(aes(y = perc_positive_cumulative-40, color = "Cases")) +
scale_y_continuous(sec.axis = sec_axis(~.*1+40, name = "Cumulative percent of tests that are positive")) +
scale_colour_manual(values = c("red", "blue")) +
labs(y = "Change in percent of devices leaving home 6 days ago relative to before", x = "Date", color = "Data", title = "San Mateo County Positive Tests and Social Distancing, 6 Day Lag")
smc_test_small %>% ggplot(
aes(x = date)) +
geom_line(aes(y = leaving6_mov, color="Leaving home")) +
geom_line(aes(y = perc_positive_cumulative_mov-40, color = "Cases")) +
scale_y_continuous(sec.axis = sec_axis(~.*1+40, name = "Cumulative percent of tests that are positive, 7 day moving average")) +
scale_colour_manual(values = c("red", "blue")) +
labs(y = "Change in percent of devices leaving home 6 days ago relative to before, 7 day moving average", x = "Date", color = "Data", title = "San Mateo County Positive Tests and Social Distancing, 6 Day Lag")
smc_test_small %>% ggplot(
aes(x = date)) +
geom_line(aes(y = leaving6_mov, color="Leaving home")) +
geom_line(aes(y = perc_positive_movavg*100-30, color = "Cases")) +
scale_y_continuous(sec.axis = sec_axis(~.*1-30/100, name = "Percent of tests that are positive, 7 day moving average")) +
scale_colour_manual(values = c("red", "blue")) +
labs(y = "Change in percent of devices leaving home 6 days ago relative to before, 7 day moving average", x = "Date", color = "Data", title = "San Mateo County Positive Tests and Social Distancing, 6 Day Lag")
smc_nolag <- smc_test_small %>%
dplyr::select(cases_growth_daily_mov, leaving_mov) %>%
cbind(`Lag` = "0")
smc_4lag <- smc_test_small %>%
dplyr::select(cases_growth_daily_mov, leaving4_mov) %>%
cbind(`Lag` = "4")
colnames(smc_4lag)[2] = "leaving_mov"
smc_5lag <- smc_test_small %>%
dplyr::select(cases_growth_daily_mov, leaving5_mov) %>%
cbind(`Lag` = "5")
colnames(smc_5lag)[2] = "leaving_mov"
smc_6lag <- smc_test_small %>%
dplyr::select(cases_growth_daily_mov, leaving6_mov) %>%
cbind(`Lag` = "6")
colnames(smc_6lag)[2] = "leaving_mov"
smc_7lag <- smc_test_small %>%
dplyr::select(cases_growth_daily_mov, leaving7_mov) %>%
cbind(`Lag` = "7")
colnames(smc_7lag)[2] = "leaving_mov"
smc_8lag <- smc_test_small %>%
dplyr::select(cases_growth_daily_mov, leaving8_mov) %>%
cbind(`Lag` = "8")
colnames(smc_8lag)[2] = "leaving_mov"
smc_9lag <- smc_test_small %>%
dplyr::select(cases_growth_daily_mov, leaving9_mov) %>%
cbind(`Lag` = "9")
colnames(smc_9lag)[2] = "leaving_mov"
smc_10lag <- smc_test_small %>%
dplyr::select(cases_growth_daily_mov, leaving10_mov) %>%
cbind(`Lag` = "10")
colnames(smc_10lag)[2] = "leaving_mov"
smc_scatter <- rbind(smc_nolag, smc_4lag, smc_5lag, smc_6lag, smc_7lag, smc_8lag, smc_9lag, smc_10lag)
fig_cases_growth_smc <-
plot_ly (smc_scatter) %>%
add_trace(
x = ~leaving_mov,
y = ~cases_growth_daily_mov,
frame = ~`Lag`,
type = 'scatter',
mode = 'markers',
showlegend = F
) %>%
animation_button(visible = F) %>%
layout(xaxis = list(title = '% difference in leaving home'), yaxis = list(title = 'daily case growth rate'), title = list(title = "SMC"))
fig_cases_growth_smc